function h = jcb_plot_peak_positions(filename,spectrum,savefigure,fileformat)
%% jcb_plot_peak_positions.m
% a script to extract the peak wavelengths from the files output by
% lumerical fdtd, then plot their positions as a function of crystal
% position

extra_title_string='';
if nargin<4
    fileformat='pdf';
end
if nargin<3
    savefigure=false;
end
if nargin<2
    spectrum='all';
end
if nargin<1 || isempty(filename)
    % get user to pick files
    [filelist,pathname] = jcb_pick_files('*data.txt', ...
    'Please select files for which to plot peaks',true,pwd);
else
    if strcmpi(filename,'all')
        filename='b';
    end
    [filelist,pathname] = jcb_select_files(filename,'spectra');
end

filelist=sort(filelist)'; % gets them into ascii order, so locations run 1-6
filelist=filelist(:);     % column cell-vector
%%
if ~isempty(filelist)    
    n = size(filelist,1);
    peak=zeros(n,3);    % peaks, in format (scattering, absorption, extinction)
    x=zeros(n,1);       % positions to plot points on the graph
    locstr=cell(n,1);
    num_no_cryst = 0;
    num_corner = 0;
    %% get peaks
    for ii=1:n
        filename=filelist{ii};
        disp(['extracting peaks from file ' filename]);
        data=load(fullfile(pathname,filename));
        % columns 2 to 4 are scattering, absorption & extinction
        [val ind] = max(data(:,2:4));
        peak(ii,:)=data(ind,1)';

        fileinfo=regexp(filename,'dose(?<btl>\d+)_(?<btw>\d+)_(?<gap>\d+)','names');
        loc=regexp(filename,'loc(?<loc>\d+)sweep_d(?<d>\d+)','names');
        if isempty(loc)
            loc=regexp(filename,'loc(?<loc>\d+)','names');
        end
        switch str2double(loc.loc)
            case 1
                if length(fields(loc))>1
                    locstr{ii}=loc.d;
                    x(ii)=str2double(loc.d);
                else
                    locstr{ii}='10';
                    x(ii)=10;
                end
            case 2
                locstr{ii}='20';
                x(ii)=20;
            case 3
                locstr{ii}='40';
                x(ii)=40;
            case 4
                locstr{ii}='108';
                x(ii)=108;
            case 5
                if length(fields(loc))>1
                    locstr{ii}=['corner, ' loc.d];
                else
                    locstr{ii}='corner';
                end
                x(ii)=25;
                num_corner = num_corner+1;
            case 6
                locstr{ii}='no crystal';  % as no crystal used
                x(ii)=120;
                num_no_cryst = num_no_cryst+1;
        end
    end
    
    %% plot 
    h = figure(63456);
    clf;
    hold on
    xlabel('crystal position');
    ylabel('peak wavelength');
    xlim([10 125]);
    ylim([640 720]);
    legendstr={};
    if (strcmp(spectrum,'scattering') || strcmp(spectrum,'all'))
        plot(x(1:(n-num_corner-num_no_cryst)),peak(1:(n-num_corner-num_no_cryst),1),'r.-');
        plot(x(n-num_corner-num_no_cryst+1:n-num_no_cryst),peak(n-num_corner-num_no_cryst+1:n-num_no_cryst,1),'rd');
        plot(x(n-num_no_cryst+1:n),peak(n-num_no_cryst+1:n,1),'rx');
        if n-num_corner-num_no_cryst>0
            legendstr{end+1}='scattering';
        end
        if num_corner
            legendstr{end+1}='corner scattering';
        end
        if num_no_cryst
            legendstr{end+1}='no-crystal scattering';
        end
    end
    if (strcmp(spectrum,'absorption') || strcmp(spectrum,'all'))
        plot(x(1:(n-num_corner-num_no_cryst)),peak(1:(n-num_corner-num_no_cryst),2),'g.-');
        plot(x(n-num_corner-num_no_cryst+1:n-num_no_cryst),peak(n-num_corner-num_no_cryst+1:n-num_no_cryst,2),'gd');
        plot(x(n-num_no_cryst+1:n),peak(n-num_no_cryst+1:n,2),'gx');
        if n-num_corner-num_no_cryst>0
            legendstr{end+1}='absorption';
        end
        if num_corner
            legendstr{end+1}='corner absorption';
        end
        if num_no_cryst
            legendstr{end+1}='no-crystal absorption';
        end
    end
    if (strcmp(spectrum,'extinction') || strcmp(spectrum,'all'))
        plot(x(1:(n-num_corner-num_no_cryst)),peak(1:(n-num_corner-num_no_cryst),3),'b.-');
        plot(x(n-num_corner-num_no_cryst+1:n-num_no_cryst),peak(n-num_corner-num_no_cryst+1:n-num_no_cryst,3),'bd');
        plot(x(n-num_no_cryst+1:n),peak(n-num_no_cryst+1:n,3),'bx');
        if n-num_corner-num_no_cryst>0
            legendstr{end+1}='extinction';
        end
        if num_corner
            legendstr{end+1}='corner extinction';
        end
        if num_no_cryst
            legendstr{end+1}='no-crystal extinction';
        end
    end
        
    legend(legendstr);
    title(['spectra for bowties of size ' fileinfo.btl 'nm' extra_title_string]);
    hold off
	
    if savefigure
        print(['-d' fileformat],fullfile(pathname,['dose_' fileinfo.btl '_spectral_peaks']));
    end
    
end

if nargout<1
    clear h
end
end